ON THE STABILITY OF THE BAREISS AND RELATED TOEPLITZ FACTORIZATION ALGORITHMS^{∗}

A. W. BOJANCZYK^{†}, R. P. BRENT^{‡}, F. R. DE HOOG^{§}, AND D. R. SWEET ^{¶}

Abstract. This paper contains a numerical stability analysis of factorization algorithms for computing the Cholesky decomposition of symmetric positive definite matrices of displacement rank 2. The algorithms in the class can be expressed as sequences of elementary downdating steps. The stability of the factorization algorithms follows directly from the numerical properties of algorithms for realizing elementary downdating operations. It is shown that the Bareiss algorithm for factorizing a symmetric positive definite Toeplitz matrix is in the class and hence the Bareiss algorithm is stable. Some numerical experiments that compare behavior of the Bareiss algorithm and the Levinson algorithm are presented. These experiments indicate that in general (when the reflection coefficients are not all of the same sign) the Levinson algorithm can give much larger residuals than the Bareiss algorithm.

Key words. Toeplitz matrices, Bareiss algorithm, Levinson algorithm, numerical stability

AMS subject classifications. 65F05, 65G05, 47B35, 65F30

1. Introduction. We consider the numerical stability of algorithms for solving a linear system

T x = b, (1.1)

where T is an n × n positive definite Toeplitz matrix and b is an n × 1 vector. We assume that the system is solved in floating point arithmetic with relative precision ϵ by first computing the Cholesky factor of T. Hence the emphasis of the paper is on factorization algorithms for the matrix T.

Roundoff error analyses of Toeplitz systems solvers have been given by Cybenko [10] and Sweet [22]. Cybenko showed that the Levinson-Durbin algorithm produces a residual which, under the condition that all reflection coefficients are positive, is of comparable size to that produced by the well behaved Cholesky method. He hypothesised that the same is true even if the reflection coefficients are not all positive. If correct, this would indicate that numerical quality of the Levinson-Durbin algorithm is comparable to that of the Cholesky method.

In his PhD thesis [22], Sweet presented a roundoff error analysis of a variant of the Bareiss algorithm [2], and concluded that the algorithm is numerically stable (in the sense specified in Section 7). In this paper we strengthen and generalize these early results on the stability of the Bareiss algorithm. In particular, our approach via elementary downdating greatly simplifies roundoff error analysis and makes it applicable to a larger-than-Toeplitz class of matrices.

After introducing the notation and the concept of elementary downdating in Sections 2 and 3, in Section 4 we derive matrix factorization algorithms as a sequence

∗ To appear in SIAM J. Matrix Anal. Appl. (accepted Sept. 1993).

† School of Electrical Engineering, Cornell University, Ithaca, NY 14853-5401, USA,

adamb@toeplitz.ee.cornell.edu

‡ Computer Sciences Laboratory, Australian National University, Canberra, ACT 0200, Australia,

rpb@cslab.anu.edu.au

§ Division of Mathematics and Statistics, CSIRO, Canberra, ACT 2601, Australia,

frank@dmscanb.cbr.dms.csiro.au

¶ Electronics Research Lab., DSTO, Salisbury, SA 5108, Australia, dougs@ewd.dsto.gov.au Copyright \begin{matrix}c \\ ◯ \end{matrix} 1993, the authors.

rpb144 typeset using LaTEX

of elementary downdating operations (see also [4]). In Section 5 we present a first order analysis by bounding the first term in an asymptotic expansion for the error in powers of ϵ. By analyzing the propagation of first order error in the sequence of downdatings that define the algorithms, we obtain bounds on the perturbations of the factors in the decompositions. We show that the computed upper triangular factor Ũ of a positive definite Toeplitz matrix T satisfies

T = Ũ^{T} Ũ + ΔT , ||ΔT || ≤ c(n)ϵ||T || ,

where c(n) is a low order polynomial in n and is independent of the condition number of T. Many of the results of Sections 2–5 were first reported in [5], which also contains some results on the stability of Levinson’s algorithm.

In Section 6 we discuss the connection with the Bareiss algorithm and conclude that the Bareiss algorithm is stable for the class of symmetric positive definite matrices. Finally, in Section 7 we report some interesting numerical examples that contrast the behaviour of the Bareiss algorithm with that of the Levinson algorithm. We show numerically that, in cases where the reflection coefficients are not all of the same sign, the Levinson algorithm can give much larger residuals than the Bareiss or Cholesky algorithms.

2. Notation. Unless it is clear from the context, all vectors are real and of dimension n. Likewise, all matrices are real and their default dimension is n × n. If a ∈ ℜ^{n}, ∥a∥ denotes the usual Euclidean norm, and if T ∈ ℜ^{n×n}, ∥T ∥ denotes the induced matrix norm:

∥T ∥ = max_{∥a∥=1} ∥T a∥ .

Our primary interest is in a symmetric positive definite Toeplitz matrix T whose i, jth entry is

t_{ij} = t_{|i−j|}.

We denote by e_{k}, k = 1, …, n, the unit vector whose kth element is 1 and whose other elements are 0. We use the following special matrices:

Z ≡ ∑_{i=1}^{n−1} e_{i+1}e_{i}^{T} = ( \begin{matrix}0 & & ⋯ & ⋯ & 0 \\ 1 & 0 & ⋯ & 0 & 0 \\ 0 & ... & & ⋮ & ⋮ \\ ⋮ & & ... & 0 & 0 \\ 0 & ⋯ & 0 & 1 & 0 \end{matrix} ) ,

J ≡ ∑_{i=1}^{n} e_{n−i+1}e_{i}^{T} = ( \begin{matrix}0 & ⋯ & ⋯ & 0 & 1 \\ ⋮ & & · & 1 & 0 \\ ⋮ & · & · & · & ⋮ \\ 0 & 1 & · & & ⋮ \\ 1 & 0 & ⋯ & ⋯ & 0 \end{matrix} ) .

The matrix Z is known as a shift-down matrix. We also make use of powers of the matrix Z, for which we introduce the following notation:

Z_{k} = \{ \begin{matrix}I & if k = 0, \\ Z^{k} & if k > 0. \end{matrix}

The antidiagonal matrix J is called a reversal matrix, because the effect of applying J to a vector is to reverse the order of components of the vector:

J [ x_{1} ; x_{2} ; ⋮ ; x_{n} ] = [ x_{n} ; x_{n−1} ; ⋮ ; x_{1} ] .

The hyperbolic rotation matrix H(θ) ∈ ℜ^{2×2} is defined by

(2.1) H(θ) = \frac{1}{cos θ} [ \begin{matrix}1 & − sin θ \\ − sin θ & 1 \end{matrix} ] .

The matrix H(θ) satisfies the relation

H(θ) [ \begin{matrix}1 & 0 \\ 0 & −1 \end{matrix} ] H(θ) = [ \begin{matrix}1 & 0 \\ 0 & −1 \end{matrix} ] ,

and it has eigenvalues λ_{1}(θ), λ_{2}(θ) given by

λ_{1}(θ) = λ_{2}^{−1}(θ) = sec θ − tan θ. (2.2)

For a given pair of real numbers a and b with |a| > |b|, there exists a hyperbolic rotation matrix H(θ) such that

H(θ) [ \begin{matrix}a \\ b \end{matrix} ] = [ \begin{matrix}√ & \overline{a^{2} − b^{2}} \\ 0 \end{matrix} ] . (2.3)

The angle of rotation θ is determined by

sin θ = \frac{b}{a} , cos θ = \frac{\sqrt{a^{2} − b^{2}}}{a} . (2.4)

3. Elementary downdating. In this section we introduce the concept of elementary downdating. The elementary downdating problem is a special case of a more general downdating problem that arises in Cholesky factorization of a positive definite difference of two outer product matrices [1, 6, 7, 12]. In Section 4, factorization algorithms are derived in terms of a sequence of downdating steps. The numerical properties of the algorithms are then related to the properties of the sequence of elementary downdating steps.

Let u_{k}, v_{k} ∈ ℜ^{n} have the following form:

k
↓
u_{k}^{T} = [0 … 0 × × × … × ] ,
v_{k}^{T} = [0 … 0 0  × × … × ] ,
↑
k + 1

that is:

e_{j}^{T} u_{k} = 0 , j < k , and e_{j}^{T} v_{k} = 0 , j ≤ k .

Applying the shift-down matrix Z to u_{k}, we have

k + 1

↓
u_{k}^{T} Z^{T} = [0 … 0 0 × × … ×] ,
v_{k}^{T}  = [0 … 0 0 × × … ×] .
↑
k + 1

Suppose that we wish to find u_{k+1}, v_{k+1} ∈ ℜ^{n} to satisfy

(3.1) u_{k+1}u_{k+1}^{T} − v_{k+1}v_{k+1}^{T} = Zu_{k}u_{k}^{T} Z^{T} − v_{k}v_{k}^{T} ,

where

k + 1
↓
u_{k+1}^{T} = [0 … 0 0 ×  × … ×] ,
v_{k+1}^{T} = [0 … 0 0 0 × … ×] ,
↑
k + 2

that is

e_{j}^{T} u_{k+1} = 0 , j < k + 1 , and e_{j}^{T} v_{k+1} = 0 , j ≤ k + 1 . (3.2)

We refer to the problem of finding u_{k+1} and v_{k+1} to satisfy (3.1), given u_{k} and v_{k}, as the elementary downdating problem. It can be rewritten as follows:

[u_{k+1} v_{k+1}] [ \begin{matrix}1 & 0 \\ 0 & −1 \end{matrix} ] [ \begin{matrix}u_{k+1}^{T} \\ v_{k+1}^{T} \end{matrix} ] = [Zu_{k} v_{k}] [ \begin{matrix}1 & 0 \\ 0 & −1 \end{matrix} ] [ \begin{matrix}u_{k}^{T} Z^{T} \\ v_{k}^{T} \end{matrix} ] .

From (2.1), (2.3) and (2.4), it is clear that the vectors u_{k+1} and v_{k+1} can be found by using a hyperbolic rotation H (θ_{k}) defined by the following relations:

(3.3a)  sin θ_{k} = e_{k+1}^{T}v_{k}/e_{k}^{T} u_{k} ,
(3.3b)  cos θ_{k} = \sqrt{1 − sin^{2} θ_{k}} ,

and

[ \begin{matrix}u_{k+1}^{T} \\ v_{k+1}^{T} \end{matrix} ] = H (θ_{k}) [ \begin{matrix}u_{k}^{T} Z^{T} \\ v_{k}^{T} \end{matrix} ] . (3.4)

The elementary downdating problem has a unique solution (up to sign changes) if

|e_{k}^{T} u_{k}| > |e_{k+1}^{T}v_{k}| .

The calculation of u_{k+1}, v_{k+1} via (3.4) can be performed in the obvious manner. Following common usage, algorithms which perform downdating in this manner will be referred to as hyperbolic downdating algorithms.

Some computational advantages may be obtained by rewriting (3.1) as follows:

[u_{k+1} v_{k}] [ \begin{matrix}u_{k+1}^{T} \\ v_{k}^{T} \end{matrix} ] = [Zu_{k} v_{k+1}] [ \begin{matrix}u_{k}^{T} Z^{T} \\ v_{k+1}^{T} \end{matrix} ] .

Consider now an orthogonal rotation matrix G(θ_{k}),

G(θ_{k}) = [ \begin{matrix}cos θ_{k} & sin θ_{k} \\ − sin θ_{k} & cos θ_{k} \end{matrix} ] ,

where cos θ_{k} and sin θ_{k} are defined by (3.3b) and (3.3a), respectively. Then it is easy to check that

G(θ_{k}) [ \begin{matrix}u_{k+1}^{T} \\ v_{k}^{T} \end{matrix} ] = [ \begin{matrix}u_{k}^{T} Z^{T} \\ v_{k+1}^{T} \end{matrix} ] , (3.5)

or, equivalently,

[ \begin{matrix}u_{k+1}^{T} \\ v_{k}^{T} \end{matrix} ] = G(θ_{k})^{T} [ \begin{matrix}u_{k}^{T} Z^{T} \\ v_{k+1}^{T} \end{matrix} ] . (3.6)

Thus, we may rewrite (3.6) as

(3.7a)  v_{k+1} = (v_{k} − sin θ_{k}Zu_{k})/ cos θ_{k} ,
(3.7b)  u_{k+1} = − sin θ_{k}v_{k+1} + cos θ_{k}Zu_{k} .

Note that equation (3.7a) is the same as the second component of (3.4). However, (3.7b) differs from the first component of (3.4) as it uses v_{k+1} in place of v_{k} to define u_{k+1}. It is possible to construct an alternative algorithm by using the first component of (3.5) to define u_{k+1}. This leads to the following formulas:

(3.8a)  u_{k+1} = (Zu_{k} − sin θ_{k}v_{k})/ cos θ_{k} ,
(3.8b)  v_{k+1} = − sin θ_{k}u_{k+1} + cos θ_{k}v_{k} .

We call algorithms based on (3.7a)–(3.7b) or (3.8a)–(3.8b) mixed elementary downdating algorithms. The reason for considering mixed algorithms is that they have superior stability properties to hyperbolic algorithms in the following sense.

Let ũ_{k}, ṽ_{k} be the values of u_{k}, v_{k} that are computed in floating point arithmetic with relative machine precision ϵ. The computed values ũ_{k}, ṽ_{k} satisfy a perturbed version of (3.1), that is,

(3.9)  ũ_{k+1}ũ_{k+1}^{T} − ṽ_{k+1}ṽ_{k+1}^{T} = Zũ_{k}ũ_{k}^{T} Z^{T} − ṽ_{k}ṽ_{k}^{T} + ϵG_{k} + O(ϵ^{2}) ,

where the second order term O(ϵ^{2}) should be understood as a matrix whose elements are bounded by a constant multiple of ϵ^{2}||G_{k}||. The norm of the perturbation G_{k} depends on the precise specification of the algorithm used. It can be shown [6] that the term G_{k} satisfies

∥G_{k}∥ ≤ c_{m} (∥Zu_{k}∥^{2} + ∥v_{k}∥^{2} + ∥u_{k+1}∥^{2} + ∥v_{k+1}∥^{2}) (3.10)

when a mixed downdating strategy is used (here c_{m} is a positive constant). When hyperbolic downdating is used the term G_{k} satisfies

∥G_{k}∥ ≤ c_{h}∥H(θ_{k})∥ (∥Zu_{k}∥ + ∥v_{k}∥) (∥u_{k+1}∥ + ∥v_{k+1}∥) , (3.11)

where c_{h} is a positive constant [6]. (The constants c_{m} and c_{h} are dependent on implementation details, but are of order unity and independent of n.) Note the presence of the multiplier ∥H(θ_{k})∥ in the bound (3.11) but not in (3.10). In view of (2.2), ∥H(θ_{k})∥ could be large. The significance of the multiplier ∥H(θ_{k})∥ depends on the context in which the downdating arises. We consider the implications of the bounds (3.10) and (3.11) in Section 5 after we make a connection between downdating and the factorization of Toeplitz matrices.

It is easily seen that a single step of the hyperbolic or mixed downdating algorithm requires 4(n − k) + O(1) multiplications. A substantial increase in efficiency can be achieved by considering the following modified downdating problem. Given α_{k}, β_{k} ∈ ℜ and w_{k}, x_{k} ∈ ℜ^{n} that satisfy

e_{j}^{T} w_{k} = 0 , j < k and e_{j}^{T} x_{k} = 0 , j ≤ k ,

find α_{k+1}, β_{k+1} and w_{k+1}, x_{k+1} ∈ ℜ^{n} that satisfy

α_{k+1}^{2}w_{k+1}w_{k+1}^{T} − β_{k+1}^{2}x_{k+1}x_{k+1}^{T} = α_{k}^{2}Zw_{k}w_{k}^{T} Z^{T} − β_{k}^{2}x_{k}x_{k}^{T} ,

with

e_{j}^{T} w_{k} = 0 , j < k and e_{j}^{T} x_{k} = 0 , j ≤ k .

If we make the identification

u_{k} = α_{k}w_{k} and v_{k} = β_{k}x_{k} ,

then we find that the modified elementary downdating problem is equivalent to the elementary downdating problem. However, the extra parameters can be chosen judiciously to eliminate some multiplications. For example, if we take α_{k} = β_{k}, α_{k+1} = β_{k+1}, then from (3.3a), (3.3b) and (3.4),

(3.12a)  sin θ_{k} = e_{k+1}^{T}x_{k}/e_{k}^{T} w_{k} ,
(3.12b)  α_{k+1} = α_{k}/ cos θ_{k} ,

and

(3.13a)  w_{k+1} = Zw_{k} − sin θ_{k}x_{k} ,
(3.13b)  x_{k+1} = − sin θ_{k}Zw_{k} + x_{k} .

Equations (3.12a)–(3.13b) form a basis for a scaled hyperbolic elementary downdating algorithm which requires 2(n − k) + O(1) multiplications. This is about half the number required by the unscaled algorithm based on (3.4). (The price is an increased likelihood of underflow or overflow, but this can be avoided if suitable precautions are taken in the code.)

Similarly, from (3.7a) and (3.7b) we can obtain a scaled mixed elementary downdating algorithm via

sin θ_{k} = β_{k}e_{k+1}^{T}x_{k}/α_{k}e_{k}^{T} w_{k} ,
α_{k+1} = α_{k} cos θ_{k} ,
β_{k+1} = β_{k}/ cos θ_{k} ,

and

x_{k+1} = x_{k} − \frac{sin θ_{k}α_{k}}{β_{k}} Zw_{k} ,
w_{k+1} = −\frac{sin θ_{k}β_{k+1}}{α_{k+1}} x_{k+1} + Zw_{k} .

The stability properties of scaled mixed algorithms are similar to those of the corresponding unscaled algorithms [12].

4. Symmetric factorization. We adopt the following definition from [18]. Definition 4.1. An n × n symmetric matrix T has displacement rank 2 iff there

exist vectors u, v ∈ ℜ^{n} such that

T − ZT Z^{T} = uu^{T} − vv^{T} . (4.1)

The vectors u and v are called the generators of T and determine the matrix T uniquely. Whenever we want to stress the dependence of T on u and v we write T = T (u , v).

In the sequel we will be concerned with a subset T of all matrices satisfying (4.1). The subset is defined as follows.

Definition 4.2. A matrix T is in T iff

(a) T is positive definite,

(b) T satisfies (4.1) with generators u and v,

(c) v^{T} e_{1} = 0, i.e., the first component of v is zero.

It is well known that positive definite n × n Toeplitz matrices form a subset of T . Indeed, if T = (t_{|i−j|})_{i,j=0}^{n−1}, then

T − ZT Z^{T} = uu^{T} − vv^{T} ,

where

u^{T} = (t_{0} , t_{1} , …, t_{n−1}) /\sqrt{t_{0}} ,
v^{T} = (0 , t_{1} , …, t_{n−1}) /\sqrt{t_{0}} .

The set T also contains matrices which are not Toeplitz, as the following example shows. Example. Let

T = [ \begin{matrix}25 & 20 & 15 \\ 20 & 32 & 29 \\ 15 & 29 & 40 \end{matrix} ] , u = [ \begin{matrix}5 \\ 4 \\ 3 \end{matrix} ] and v = [ \begin{matrix}0 \\ 3 \\ 1 \end{matrix} ] .

It is easy to check that T is positive definite. Moreover,

T − ZT Z^{T} = [ \begin{matrix}25 & 20 & 15 \\ 20 & 7 & 9 \\ 15 & 9 & 8 \end{matrix} ] = [ \begin{matrix}25 & 20 & 15 \\ 20 & 16 & 12 \\ 15 & 12 & 9 \end{matrix} ] − [ \begin{matrix}0 & 0 & 0 \\ 0 & 9 & 3 \\ 0 & 3 & 1 \end{matrix} ] = uu^{T} − vv^{T} .

Hence T = T (u , v) ∈ T , but T is not Toeplitz.

We now establish a connection between the elementary downdating problem and symmetric factorizations of a matrix from the set T .

Let T = T (u , v) ∈ T . Set

u_{1} = u, v_{1} = v

and, for k = 1, …, n − 1, solve the elementary downdating problem defined by (3.1),

u_{k+1}u_{k+1}^{T} − v_{k+1}v_{k+1}^{T} = Zu_{k}u_{k}^{T} Z^{T} − v_{k}v_{k}^{T} ,

which we assume for the moment has a solution for each k. On summing over k = 1, …, n − 1 we obtain

∑_{k=1}^{n−1} u_{k+1}u_{k+1}^{T} − ∑_{k=1}^{n−1} v_{k+1}v_{k+1}^{T} = ∑_{k=1}^{n−1} Zu_{k}u_{k}^{T} Z^{T} − ∑_{k=1}^{n−1} v_{k}v_{k}^{T} .

If we now observe that, from (3.2),

Zu_{n} = v_{n} = 0 ,

we arrive at the following relation:

(4.2)  ∑_{k=1}^{n} u_{k}u_{k}^{T} − Z (∑_{k=1}^{n} u_{k}u_{k}^{T} ) Z^{T} = u_{1}u_{1}^{T} − v_{1}v_{1}^{T} ,

which implies that ∑_{k=1}^{n} u_{k}u_{k}^{T} ∈ T . Moreover, as matrices having the same generators are identical, we obtain

T = ∑_{k=1}^{n} u_{k}u_{k}^{T} = U^{T} U ,

where

U = ∑_{k=1}^{n} e_{k}u_{k}^{T}

is upper triangular, and hence is the Cholesky factor of T. We have derived, albeit in a rather indirect manner, the basis of an algorithm for calculating the Cholesky decomposition of a matrix from the set T .

We now return to the question of existence of a solution to the elementary downdating problem for each k = 1, …, n − 1. It is easy to verify that, if T ∈ T , then |e_{1}^{T} u_{1}| > |e_{2}^{T} v_{1}|. Using (4.2) and (3.1), it can be shown by induction on k that

|e_{k}^{T} u_{k}| > |e_{k+1}^{T}v_{k}|, k = 2, …, n − 1.

Consequently, | sin θ_{k}| < 1 in (3.3a), and the elementary downdating problem has a solution for each k = 1, …, n − 1.

To summarize, we have the following algorithm for factorizing a matrix T = T (u , v) ∈ T .

Algorithm FACTOR( T):

Set u_{1} = u, v_{1} = v.

For k = 1, …, n − 1 calculate u_{k+1}, v_{k+1} such that
u_{k+1}u_{k+1}^{T} − v_{k+1}v_{k+1}^{T} = Zu_{k}u_{k}^{T} Z^{T} − v_{k}v_{k}^{T} ,
e_{k+1}^{T}v_{k+1} = 0 .
Then T = U^{T} U, where U = ∑_{k=1}^{n} e_{k}u_{k}^{T}.

In fact we have not one algorithm but a class of factorization algorithms, where each algorithm corresponds to a particular way of realizing the elementary downdating steps. For example, the connection with the scaled elementary downdating problem is straightforward. On making the identification

u_{k} = α_{k}w_{k} and v_{k} = β_{k}x_{k} , (4.3)

we obtain

T = W^{T} D^{2}W ,

where

W = ∑_{k=1}^{n} e_{k}w_{k}^{T} ,

D = ∑_{k=1}^{n} α_{k}e_{k}e_{k}^{T} .

It is clear from Section 3 that Algorithm FACTOR( T) requires 2n^{2} +O(n) multiplications when the unscaled version of elementary downdating is used, and n^{2} +O(n) multiplications when the scaled version of elementary downdating is used. However, in the sequel we do not dwell on the precise details of algorithms. Using (4.3), we can relate algorithms based on the scaled elementary downdating problem to those based on the unscaled elementary downdating problem. Thus, for simplicity, we consider only the unscaled elementary downdating algorithms.

5. Analysis of factorization algorithms. In this section we present a numerical stability analysis of the factorization of T ∈ T via Algorithm FACTOR( T). The result of the analysis is applied to the case when the matrix T is Toeplitz.

Let ũ_{k}, ṽ_{k} be the values of u_{k}, v_{k} that are computed in floating point arithmetic with relative machine relative precision ϵ. The computed quantities ũ_{k} and ṽ_{k} satisfy the relations

ũ_{k} = u_{k} + O(ϵ), ṽ_{k} = v_{k} + O(ϵ), (5.1)

and the aim of this section is to provide a first order analysis of the error. By a first order analysis we mean that the error can be bounded by a function which has an asymptotic expansion in powers of ϵ, but we only consider the first term of this asymptotic expansion. One should think of ϵ → 0+ while the problem remains fixed [19]. Thus, in this section (except for Corollary 5.5) we omit functions of n from the “O” terms in relations such as (5.1) and (5.2).

The computed vectors ũ_{k}, ṽ_{k} satisfy a perturbed version (3.9) of (3.1). On summing (3.9) over k = 1, …, n − 1 we obtain

T̃ − ZT̃Z^{T} = ũ_{1}ũ_{1}^{T} − ṽ_{1}ṽ_{1}^{T} − (Zũ_{n}ũ_{n}^{T} Z^{T} − ṽ_{n}ṽ_{n}^{T}) + ϵ ∑_{k=1}^{n−1} G_{k} + O(ϵ^{2}) ,

where

T̃ = Ũ^{T} Ũ ,
Ũ = ∑_{k=1}^{n} e_{k}ũ_{k}^{T} .

Since

Zũ_{n} = O(ϵ), ṽ_{n} = O(ϵ) ,

we find that

T̃ − ZT̃Z^{T} = ũ_{1}ũ_{1}^{T} − ṽ_{1}ṽ_{1}^{T} + ϵ ∑_{k=1}^{n−1} G_{k} + O(ϵ^{2}) . (5.2)

Now define

Ẽ = T̃ − T. (5.3)

Then, using (4.1), (5.2) and (5.3),

Ẽ − ZẼZ^{T} = ũ_{1}ũ_{1}^{T} − uu^{T} + ṽ_{1}ṽ_{1}^{T} − vv^{T} + ϵ ∑_{k=1}^{n−1} G_{k} + O(ϵ^{2}) .

In a similar manner we obtain expressions for Z_{j}ẼZ_{j}^{T} − Z_{j+1}ẼZ_{j+1}^{T}, j = 0, …, n − 1. Summing over j gives

Ẽ =∑_{j=0}^{n−1} Z_{j}((ũ_{1}ũ_{1}^{T} − u_{1}u_{1}^{T}) + (ṽ_{1}ṽ_{1}^{T} − v_{1}v_{1}^{T}))Z_{j}^{T}

(5.4)

+ ϵ ∑_{j=0}^{n−1} ∑_{k=1}^{n−1} Z_{j}G_{k}Z_{j}^{T} + O(ϵ^{2}) .

We see from (5.4) that the error consists of two parts – the first associated with initial errors and the second associated with the fact that (5.2) contains an inhomogeneous term. Now

∥ũ_{1}ũ_{1}^{T} − uu^{T} ∥ ≤ 2∥u∥ ∥ũ_{1} − u∥ + O(ϵ^{2}) ,
∥ṽ_{1}ṽ_{1}^{T} − vv^{T} ∥ ≤ 2∥v∥ ∥ṽ_{1} − v∥ + O(ϵ^{2}) .

Furthermore, from (4.1),

T r(T) − T r(ZT Z^{T}) = ∥u∥^{2} − ∥v∥^{2} > 0 ,

and hence

∥∑_{j=0}^{n−1} Z_{j}(ũ_{1}ũ_{1}^{T} − uu^{T} + ṽ_{1}ṽ_{1}^{T} − vv^{T})Z_{j}^{T} ∥

(5.5)

≤ 2n∥u∥(∥ũ_{1} − u∥ + ∥ṽ_{1} − v∥) + O(ϵ^{2}) .

This demonstrates that initial errors do not propagate unduly. To investigate the double sum in (5.4) we require a preliminary result.

Lemma 5.1. For k = 1, 2, …, n − 1 and j = 0, 1, 2, …,

∥Z_{j}v_{k}∥ ≤ ∥Z_{j+1}u_{k}∥ .

Proof. Let

T_{k} = T − ∑_{l=1}^{k} u_{l}u_{l}^{T} = ∑_{l=k+1}^{n} u_{l}u_{l}^{T} .

It is easy to verify that

T_{k} − ZT_{k}Z^{T} = Zu_{k}u_{k}^{T} Z^{T} − v_{k}v_{k}^{T}

and, since T_{k} is positive semi-definite,

T r(Z_{j}T_{k}Z_{j}^{T} − Z_{j+1}T_{k}Z_{j+1}^{T}) = ∥Z_{j+1}u_{k}∥^{2} − ∥Z_{j}v_{k}∥^{2} ≥ 0 .

We now demonstrate stability when the mixed version of elementary downdating is used in Algorithm FACTOR( T). In this case the inhomogeneous term G_{k} satisfies a shifted version of (3.10), that is

∥Z_{j}G_{k}Z_{j}^{T} ∥ ≤ c_{m}(∥Z_{j+1}u_{k}∥^{2} + ∥Z_{j}v_{k}∥^{2} + ∥Z_{j}u_{k+1}∥^{2} + ∥Z_{j}v_{k+1}∥^{2}) , (5.6)

where c_{m} is a positive constant.

Theorem 5.2. Assume that (3.9) and (5.6) hold. Then

∥T − Ũ^{T} Ũ∥ ≤ 2n∥u∥(∥ũ_{1} − u∥ + ∥ṽ_{1} − v∥) + 4ϵc_{m} ∑_{j=0}^{n−1} T r(Z_{j}T Z_{j}^{T}) + O(ϵ^{2}) .

Proof. Using Lemma 5.1,

∥Z_{j}G_{k}Z_{j}^{T} ∥ ≤ 2c_{m}(∥Z_{j+1}u_{k}∥^{2} + ∥Z_{j}u_{k+1}∥^{2}).

Furthermore, since

T r(Z_{j}T Z_{j}^{T}) = ∑_{k=1}^{n} ∥Z_{j}u_{k}∥^{2},

it follows that

∥∑_{j=0}^{n−1} ∑_{k=1}^{n} Z_{j}G_{k}Z_{j}^{T} ∥ ≤ 4c_{m} ∑_{j=0}^{n−1} T r(Z_{j}T Z_{j}^{T}) . (5.7)

The result now follows from (5.4), (5.5) and (5.7).

For the hyperbolic version of the elementary downdating algorithms a shifted version of the weaker bound (3.11) on G_{k} holds (see [6]), namely

∥Z_{j}G_{k}Z_{j}^{T} ∥ ≤ c_{h}∥H(θ_{k})∥(∥Z_{j+1}u_{k}∥ + ∥Z_{j}v_{k}∥)(∥Z_{j}u_{k+1}∥ + ∥Z_{j}v_{k+1}∥) . (5.8)

By Lemma 5.1, this simplifies to

∥Z_{j}G_{k}Z_{j}^{T} ∥ ≤ 4c_{h}∥H(θ_{k})∥ ∥Z_{j+1}u_{k}∥ ∥Z_{j}u_{k+1}∥ . (5.9)

The essential difference between (3.10) and (3.11) is the occurence of the multiplier ∥H(θ_{k})∥ which can be quite large. This term explains numerical difficulties in applications such as the downdating of a Cholesky decomposition [6]. However, because of the special structure of the matrix T, it is of lesser importance here, in view of the following result.

Lemma 5.3. For k = 1, 2, …, n − 1 and j = 0, 1, …, n − k,

∥H(θ_{k})∥ ∥Z_{j}u_{k+1}∥ ≤ 2(n − k − j)∥Z_{j+1}u_{k}∥.

Proof. It is easy to verify from (3.4) that

\frac{1 ∓ sin θ_{k}}{cos θ_{k}} (u_{k+1} ∓ v_{k+1}) = Zu_{k} ∓ v_{k} ,

and from (2.1) that

∥H(θ_{k})∥ = \frac{1 + | sin θ|}{cos θ} .

Thus,

∥H(θ_{k})∥ ∥Z_{j}u_{k+1}∥ ≤ ∥H(θ_{k})∥ ∥Z_{j}v_{k+1}∥ + ∥Z_{j+1}u_{k}∥ + ∥Z_{j}v_{k}∥
≤ ∥H(θ_{k})∥ ∥Z_{j+1}u_{k+1}∥ + 2∥Z_{j+1}u_{k}∥ ,

where the last inequality was obtained using Lemma 5.1. Thus

∥H(θ_{k})∥ ∥Z_{j}u_{k+1}∥ ≤ 2 ∑_{l=j+1}^{n−k} ∥Z_{l}u_{k}∥ ,

and the result follows.

Remark. Lemma 5.3 does not hold for the computed quantities unless we introduce an O(ϵ) term. However in a first order analysis we only need it to hold for the exact quantities.

Theorem 5.4. Assume that (3.9) and (5.8) hold. Then

∥T − Ũ^{T} Ũ∥ ≤ 2n∥u∥(∥ũ_{1} − u∥ + ∥ṽ_{1} − v∥) + 8ϵc_{h} ∑_{j=1}^{n−1}(n − j)T r(Z_{j}T Z_{j}^{T}) + O(ϵ^{2}) .

Proof. Applying Lemma 5.3 to (5.9) gives

∥Z_{j}G_{k}Z_{j}^{T} ∥ ≤ 8c_{h}(n − j − 1)∥Z_{j+1}u_{k}∥^{2} ,

and hence

∥∑_{j=0}^{n−1} ∑_{k=1}^{n−1} Z_{j}G_{k}Z_{j}^{T} ∥ ≤ 8c_{h} ∑_{j=1}^{n−1} ∑_{k=1}^{n−1}(n − j)∥Z_{j}u_{k}∥^{2}

(5.10)

≤ 8c_{h} ∑_{j=1}^{n−1}(n − j)T r(Z_{j}T Z_{j}^{T}) .

The result now follows from (5.4), (5.5) and (5.10).

Note that, when T is Toeplitz,

T r(Z_{j}T Z_{j}^{T}) = (n − j)t_{0} .

Hence, from Theorems 5.2 and 5.4, we obtain our main result on the stability of the factorization algorithms based on Algorithm FACTOR( T) for a symmetric positive definite Toeplitz matrix:

Corollary 5.5. The factorization algorithm FACTOR( T) applied to a symmetric positive definite Toeplitz matrix T produces an upper triangular matrix Ũ such that

T = Ũ^{T} Ũ + ΔT ,

where

∥ΔT ∥ = O(ϵt_{0}n^{2})

when mixed downdating is used, and

∥ΔT ∥ = O(ϵt_{0}n^{3})

when hyperbolic downdating is used.

6. The connection with the Bareiss algorithm. In 1969, Bareiss [2] proposed an O(n^{2}) algorithm for solving Toeplitz linear systems. For a symmetric Toeplitz matrix T, the algorithm, called a symmetric Bareiss algorithm in [22], can be expressed as follows. Start with a matrix A^{(0)} := T and partition it in two ways:

A^{(0)} = ( \begin{matrix}U^{(0)} \\ T^{(1)} \end{matrix} ) = ( \begin{matrix}T^{(−1)} \\ L^{(0)} \end{matrix} ) ,

where U^{(0)} is the first row of T and L^{(0)} is the last row of T. Now, starting from A^{(0)}, compute successively two matrix sequences {A^{(i)}} and {A^{(−i)}}, i = 1, …, n − 1, according to the relations

A^{(i)} = A^{(i−1)} − α_{i−1}Z_{i}A^{(−i+1)} ,

(6.1)

A^{(−i)} = A^{(−i+1)} − α_{−i+1}Z_{i}^{T} A^{(i−1)} .

For 1 ≤ i ≤ n − 1, partition A^{(i)} and A^{(−i)} as follows:

A^{(i)} = ( \begin{matrix}U^{(i)} \\ T^{(i+1)} \end{matrix} ) , A^{(−i)} = ( \begin{matrix}T^{(−i−1)} \\ L^{(i)} \end{matrix} ) ,

where U^{(i)} denotes the first i + 1 rows of A^{(i)}, and L^{(i)} denotes the last i + 1 rows of A^{(−i)}. It is shown in [2] that

(a) T^{(i+1)} and T^{(−i−1)} are Toeplitz,

(b) for a proper choice of α_{i−1} and α_{−i+1}, the matrices L^{(i)} and U^{(i)} are lower and upper trapezoidal, respectively, and

(c) with the choice of α_{i−1} and α_{−i+1} as in (b), the Toeplitz matrix T^{(−i−1)} has zero elements in positions 2, …, i + 1 of its first row, while the Toeplitz matrix T^{(i+1)} has zero elements in positions n − 1, …, n − i of its last row.

Pictorially,

A^{(i)} = ( \frac{U^{(i)}}{T^{(i+1)}} ) = ( \begin{matrix}× & ⋯ & ⋯ & ⋯ & ⋯ & ⋯ & ⋯ & × \\ 0 & × & & & & & & × \\ ⋮ & ... & × & × & & & & ⋮ \\ 0 & ⋯ & 0 & × & ⋯ & ⋯ & ⋯ & × \\ × & 0 & ⋯ & 0 & × & ⋯ & ⋯ & × \\ ⋮ & ... & & & ... & ... & & ⋮ \\ ⋮ & & ... & & & ... & ... & ⋮ \\ × & ⋯ & ⋯ & × & 0 & ⋯ & 0 & × \end{matrix} )

and

A^{(−i)} = ( \frac{T^{(−i−1)}}{L^{(i)}} ) = ( \begin{matrix}× & 0 & ⋯ & 0 & × & ⋯ & ⋯ & × \\ ⋮ & ... & & & ... & ... & & ⋮ \\ ⋮ & & ... & & & ... & ... & ⋮ \\ × & ⋯ & ⋯ & × & 0 & ⋯ & 0 & × \\ × & ⋯ & ⋯ & × & × & 0 & ⋯ & 0 \\ ⋮ & & & & & ... & ... & ⋮ \\ ⋮ & & & & & & ... & 0 \\ × & × & ⋯ & ⋯ & ⋯ & ⋯ & ⋯ & × \end{matrix} ) .

After n−1 steps, the matrices A^{(n−1)} and A^{(−n+1)} are lower and upper triangular, respectively. At step i only rows i + 1, …, n of A^{(i)} and rows 1, 2, …, n − i of A^{(−i)} are modified; the remaining rows stay unchanged. Moreover, Bareiss [2] noticed that, because of the symmetry of T,

T^{(i+1)} = J_{i+1}T^{(−i−1)}J_{n} and α_{i−1} = α_{−i+1} , (6.2)

Here J_{i+1} and J_{n} are the reversal matrices of dimension (i + 1) × (i + 1) and n × n respectively.

Now, taking into account (6.2), it can be seen that the essential part of a step of the Bareiss algorithm (6.1) can be written as follows:

( \begin{matrix}t_{i+2}^{(i+1)} & t_{i+3}^{(i+1)} & … & t_{n}^{(i+1)} \\ 0 & t_{i+3}^{(−i−1)} & … & t_{n}^{(−i−1)} \end{matrix} ) =

(6.3)

( \begin{matrix}1 & −α_{i−1} \\ −α_{i−1} & 1 \end{matrix} ) ( \begin{matrix}t_{i+2}^{(i)} & t_{i+3}^{(i)} & … & t_{n}^{(i)} \\ t_{i+2}^{(−i)} & t_{i+3}^{(−i)} & … & t_{n}^{(−i)} \end{matrix} ) ,

where (t_{i+2}^{(−i)}, t_{i+2}^{(−i)}, …, t_{n}^{(−i)}) are the last (n − i − 1) components of the first row of T^{(−i)}, and (t_{i+2}^{(i)}, t_{i+3}^{(i)}, …, t_{n}^{(i)}) are the last (n − i − 1) components of the first row of T^{(i)}.

Note that (6.3) has the same form as (3.13a)–(3.13b), and hence a connection between the Bareiss algorithm and algorithm F ACT OR(T) is evident. That such a connection exists was observed by Sweet [22], and later by Delosme and Ipsen [11].

Sweet [22] related a step of the Bareiss algorithm (6.3) to a step of Bennett’s downdating procedure [3]. Next, he derived the LU factorization of a Toeplitz matrix as a sequence of Bennett’s downdating steps. Finally, he estimated the forward error in the decomposition using Fletcher and Powell’s methodology [12]. This paper generalizes and presents new derivations of the results obtained in [22].

7. Numerical examples. We adopt from [17] the following definitions of forward and backward stability.

Definition 7.1. An algorithm for solving the equation (1.1) is forward stable if the computed solution x̃ satisfies

||x − x̃|| ≤ c_{1}(n) ϵ cond(T)||x̃|| ,

where cond(T) = ||T || ||T^{−1}|| is the condition number of T, and c_{1}(n) may grow at most as fast as a polynomial in n, the dimension of the system.

Definition 7.2. An algorithm for solving the equation (1.1) is backward stable if the computed solution x̃ satisfies

||T x̃ − b|| ≤ c_{2}(n)ϵ||T || ||x̃|| ,

where c_{2}(n) may grow at most as fast as a polynomial in n, the dimension of the system.

It is known that an algorithm (for solving a system of linear equations) is backward stable iff there exists a matrix ΔT such that

(T + ΔT)x̃ = b , ||ΔT || ≤ c_{3}(n)ϵ||T || ,

where c_{3}(n) may grow at most as fast as a polynomial in n.

Note that our definitions do not require the perturbation ΔT to be Toeplitz, even if the matrix T is Toeplitz. The case that ΔT is Toeplitz is discussed in [13, 24]. The reader is referred to [9, 14, 19] for a detailed treatment of roundoff analysis for general matrix algorithms.

It is easy to see that backward stability implies forward stability, but not vice versa. This is manifested by the size of the residual vector.

Cybenko [10] showed that the L_{1} norm of the inverse of a n×n symmetric positive definite Toeplitz matrix T_{n} is bounded by

max\{\frac{1}{∏_{i=1}^{n−1} cos^{2} θ_{i}} , \frac{1}{∏_{i=1}^{n−1}(1 + sin θ_{i})}\} ≤ ∥T_{n}^{−1}∥_{1} ≤ ∏_{i=1}^{n−1} \frac{1 + | sin θ_{i}|}{1 − | sin θ_{i}|} ,

where {− sin θ_{i}}_{i=1}^{n−1} are quantities called reflection coefficients. It is not difficult to pick the reflection coefficients in such a way that the corresponding Toeplitz matrix T_{n} satisfies

cond(T_{n}) ≈ 1/ϵ .

One possible way of constructing a Toeplitz matrix with given reflection coefficients {− sin θ_{i}}_{i=1}^{n−1} is by tracing the elementary downdating steps backwards.

An example of a symmetric positive definite Toeplitz matrix that can be made poorly conditioned by suitable choice of a parameter is the prolate matrix [21, 23], defined by

t_{k} = \{ \begin{matrix}2ω & if k = 0, \\ sin(2πωk)/(πk) & otherwise, \end{matrix}

where 0 ≤ ω ≤ \frac{1}{2}. For small ω the eigenvalues of the prolate matrix cluster around 0 and 1.

We performed numerical tests in which we solved systems of Toeplitz linear equations using variants of the Bareiss and Levinson algorithms, and (for comparison) the standard Cholesky method. The relative machine precision was ϵ = 2^{−53} ≈ 10^{−16}. We varied the dimension of the system from 10 to 100, the condition number of the matrix from 1 to ϵ^{−1}, the signs of reflection coefficients, and the right hand side so the magnitude of the norm of the solution vector varied from 1 to ϵ^{−1}. In each test we monitored the errors in the decomposition, in the solution vector, and the size of the residual vector.

Table 7.1

Prolate matrix, n = 21, ω = 0.25, cond = 3.19 · 10^{14}

decomp. error soln. error resid. error
Cholesky  5.09 · 10^{−1}  7.67 · 10^{−3}  1.25 · 10^{0}
Bareiss(hyp)  3.45 · 10^{0}  1.40 · 10^{−2} 8.72 · 10^{−1}
Bareiss(mixed)  2.73 · 10^{0}  1.41 · 10^{0}  1.09 · 10^{0}
Levinson  –  5.30 · 10^{0}  4.57 · 10^{3}

Table 7.2 Reflection coefficients of the same magnitude |K| but alternating signs, |K| = 0.8956680108101296, n = 41, cond = 8.5 · 10^{15}

decomp. error soln. error resid. error
Cholesky  1.72 · 10^{−1}  6.84 · 10^{−2} 3.11 · 10^{−1}
Bareiss(hyp)  2.91 · 10^{0}  2.19 · 10^{−1} 1.15 · 10^{−1}
Bareiss(mixed)  3.63 · 10^{0}  2.48 · 10^{−1} 2.47 · 10^{−1}
Levinson  –  5.27 · 10^{−1}  1.47 · 10^{5}

Table 7.3 Reflection coefficients of the same magnitude |K| but alternating signs, |K| = 0.9795872473975045, n = 92, cond = 2.77 · 10^{15}

decomp. error soln. error resid. error
Cholesky  8.51 · 10^{−1}  3.21 · 10^{−2} 4.28 · 10^{−1}
Bareiss(hyp)  8.06 · 10^{0}  1.13 · 10^{−1} 2.28 · 10^{−1}
Bareiss(mixed)  6.71 · 10^{0}  1.16 · 10^{−1} 3.20 · 10^{−1}
Levinson  –  2.60 · 10^{−1}  1.06 · 10^{5}

Let x_{B} and x_{L} denote the solutions computed by the Bareiss and Levinson algorithms. Also, let r_{B} = T x_{B} − b and r_{L} = T x_{L} − b. Then for the Bareiss algorithms we always observed that the scaled residual

s_{B} ≡ \frac{∥r_{B}∥}{ϵ∥x_{B}∥∥T ∥}

was of order unity, as small as would be expected for a backward stable method. However, we were not able to find an example which would demonstrate the superiority of the Bareiss algorithm based on mixed downdating over the Bareiss algorithm based on hyperbolic downdating. In fact, the Bareiss algorithm based on hyperbolic downdating often gave slightly smaller errors than the Bareiss algorithm based on mixed downdating. In our experiments with Bareiss algorithms, neither the norm of the error matrix in the decomposition of T nor the residual error in the solution seemed to depend in any clear way on n, although a quadratic or cubic dependence would be expected from the worst-case error bounds of Theorems 5.2, 5.4, and Corollary 5.5.

For well conditioned systems the Bareiss and Levinson algorithms behaved similarly, and gave results comparable to results produced by a general stable method (the Cholesky method). Differences between the Bareiss and Levinson algorithms were noticeable only for very ill-conditioned systems and special right-hand side vectors.

For the Levinson algorithm, when the matrix was very ill-conditioned and the norm of the solution vector was of order unity (that is, when the norm of the solution vector did not reflect the ill-conditioning of the matrix), we often observed that the scaled residual

s_{L} ≡ \frac{∥r_{L}∥}{ϵ∥x_{L}∥∥T ∥} ,

was as large as 10^{5}. Varah [23] was the first to observe this behavior of the Levinson algorithm on the prolate matrix. Higham and Pickering [16] used a search method proposed in [15] to generate Toeplitz matrices for which the Levinson algorithm yields large residual errors. However, the search never produced s_{L} larger than 5 · 10^{5}. It plausible that s_{L} is a slowly increasing function of n and 1/ϵ.

Tables 7.1–7.3 show typical behavior of the Cholesky, Bareiss and Levinson algorithms for ill-conditioned Toeplitz systems of linear equations when the norm of the solution vectors is of order unity. The decomposition error was measured for the Cholesky and Bareiss algorithms by the quotient ||T − L · L^{T} ||/(ϵ · ||T ||), where L was the computed factor of T. The solution error was measured by the quotient ||x_{comp}−x||/||x||, where x_{comp} was the computed solution vector. Finally, the residual error was measured by the quotient ||T · x_{comp} − b||/(||T || · ||x_{comp}|| · ϵ).

8. Conclusions. This paper generalizes and presents new derivations of results obtained earlier by Sweet [22]. The bound in Corollary 5.5 for the case of mixed downdating is stronger than that given in [22]. The applicability of the Bareiss algorithms based on elementary downdating steps is extended to a class of matrices, satisfying Definition 4.2, which includes symmetric positive definite Toeplitz matrices. The approach via elementary downdating greatly simplifies roundoff error analysis. Lemmas 5.1 and 5.3 appear to be new. The stability of the Bareiss algorithms follows directly from these Lemmas and the results on the roundoff error analysis for elementary downdating steps given in [6].

The approach via downdating can be extended to the symmetric factorization of positive definite matrices of displacement rank k ≥ 2 (satisfying additional conditions similar to those listed in Definition 4.2); see [18]. For matrices of displacement rank k the factorization algorithm uses elementary rank-k downdating via hyperbolic Householder or mixed Householder reflections [8, 20].

We conclude by noting that the Bareiss algorithms guarantee small residual errors in the sense of Definition 7.2, but the Levinson algorithm can yield residuals at least five orders of magnitude larger than those expected for a backward stable method. This result suggests that, if the Levinson algorithm is used in applications where the reflection coefficients are not known in advance to be positive, the residuals should be computed to see if they are acceptably small. This can be done in O(n log n) arithmetic operations (using the FFT).

It is an interesting open question whether the Levinson algorithm can give scaled residual errors which are arbitrarily large (for matrices which are numerically nonsingular). A related question is whether the Levinson algorithm, for positive definite Toeplitz matrices T without a restriction on the reflection coefficients, is stable in the sense of Definitions 7.1 or 7.2.

REFERENCES

[1] S. T. Alexander, C-T. Pan and R. J. Plemmons, Analysis of a recursive least squares hyperbolic rotation algorithm for signal processing, Linear Algebra Appl., 98 (1988), pp. 3–40.

[2] E. H. Bareiss, Numerical solution of linear equations with Toeplitz and vector Toeplitz matrices, Numer. Math., 13 (1969), pp. 404–424.

[3] J. M. Bennett, Triangular factorization of modified matrices, Numer. Math., 7 (1965), pp. 217–221.

[4] A. W. Bojanczyk, R. P. Brent and F. R. de Hoog, QR factorization of Toeplitz matrices, Numer. Math., 49 (1986), pp. 81–94.

[5] , Stability analysis of fast Toeplitz linear system solvers, Report CMA-MR17-91, Centre for Mathematical Analysis, The Australian National University, August 1991.

[6] A. W. Bojanczyk, R. P. Brent, P. Van Dooren and F. R. de Hoog, A note on downdating the Cholesky factorization, SIAM J. Sci. Statist. Comput., 8 (1987), pp. 210–220.

[7] A. W. Bojanczyk and A. O. Steinhardt, Matrix downdating techniques for signal processing, Proc. SPIE, Vol. 975: Advanced Algorithms and Architectures for Signal Processing III, 1988, pp. 68–75.

[8] , Stabilized hyperbolic Householder transformations, IEEE Trans. Acoustics, Speech and Signal Processing, ASSP–37 (1989), pp. 1286–1288.

[9] J. R. Bunch, The weak and strong stability of algorithms in numerical linear algebra, Linear Algebra Appl., 88/89 (1987), pp. 49–66.

[10] G. Cybenko, The numerical stability of the Levinson-Durbin algorithm for Toeplitz systems of equations, SIAM J. Sci. Statist. Comput., 1 (1980), pp. 303–319.

[11] J-M. Delosme and I. C. F. Ipsen, From Bareiss’s algorithm to the stable computation of partial correlations, J. Computational and Applied Mathematics, 27 (1989), pp. 53–91.

[12] R. Fletcher and M. J. D. Powell, On the modification of LDL^{T} factorizations, Math. Comp., 28 (1974), pp. 1067–87.

[13] I. Gohberg, I. Koltracht and D. Xiao, On the solution of the Yule-Walker equation, Proc. SPIE, Vol. 1566: Advanced Algorithms and Architectures for Signal Processing IV, 1991.

[14] G. H. Golub and C. F. Van Loan, Matrix Computations, Second ed., The Johns Hopkins University Press, Baltimore, MD, 1989.

[15] N. J. Higham, Optimization by direct search in matrix computations, Numerical Analysis Report No. 197, University of Manchester, England, 1991. To appear in SIAM J. Matrix Anal. Appl.

[16] N. J. Higham and R. L. Pickering, private communication.

[17] M. Jankowski and H. Wozniakowski, Iterative refinement implies numerical stability, BIT, 17 (1977), pp. 303–311.

[18] T. Kailath, S. Y. Kung and M. Morf, Displacement ranks of matrices and linear equations, J. Math. Anal. Appl., 68 (1979), pp. 395–407.

[19] W. Miller and C. Wrathall, Software for Roundoff Analysis of Matrix Algorithms, Academic Press, New York, 1980.

[20] C. M. Rader and A. O. Steinhardt, Hyperbolic Householder transformations, IEEE Trans. Acoustics, Speech and Signal Processing, ASSP–34 (1986), pp. 1584–1602.

[21] D. Slepian, Prolate spheroidal wave functions, Fourier analysis, and uncertainty V: the discrete case, Bell Systems Tech. J., 57 (1978), pp. 1371–1430.

[22] D. R. Sweet, Numerical methods for Toeplitz matrices, Ph.D. thesis, University of Adelaide, Adelaide, Australia, 1982.

[23] J. M. Varah, The prolate matrix: a good Toeplitz test example, SIAM Conference on Control, Signals and Systems, San Francisco, 1990. Also Linear Algebra Appl., 187 (1993), pp. 269–278.

[24] , Backward error estimates for Toeplitz and Vandermonde systems, preprint, 1992. Also Tech. Report 91-20, Univ. of British Columbia, Sept. 1991.